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Abstract 

The review presents a parameter switching algorithm and his applications which allows numerical approximation of 
any attractor of a class of continuous-time dynamical systems depending linearly on a real parameter. The considered 
classes of systems are modeled by a general initial value problem embedding dynamical systems which are continuous and 
discontinuous with respect to the state variable, and of integer and fractional order. The numerous results, presented 
in several papers, are systematized here on four representative known examples representing the four classes. The 
analytical proof of the algorithm convergence for the systems belonging to the continuous class is briefly presented, 
while for the other categories of systems the convergence is numerically verified via computational tools. The utilized 
numerical tools necessary to apply the algorithm are contained in five appendices. 

Keywords Parameter switching; Global attractors; Local attractors; Fractional systems; Discontinuous systems; Filippov 
regularization 

1 Introduction 

In Nature there are many different interactions and the real systems could evolve according to more that one dynamics for 
short periods of time. Therefore, it is reasonable to think that the evolution of some natural processes could be imagined 
as the result of the alternation of different dynamics for relatively short periods of time. In particular, a topic of research 
regarding parameter switching which has arisen in the last years, consists in studying the dynamics of continuous-time 
systems [HElElllllSlEllI] and discrete systems piQllTO]. 

In this paper we review aspects of some previous results, obtained by us with parameter switching techniques. We 
considered general classes of systems, continuous and discontinuous with respect to the state variable, and of integer and 
fractional order. 

The truthfulness of the previous results are sustained here by another numerical tool, the cross-correlation. 

Via numerical simulations we found for a large class of systems (and analytically proved for a particular class [TT]). 
that any attractor of some considered system can be synthesized (numerically approximated) by using some parameter 
switching rule. The analytical and numerical proofs are based on the fact that the invariant sets obtained with the 
control parameter periodically switched are numerical approximations of those corresponding to the control parameter 
replaced with the average of the switched values. This useful result was intensively verified on several examples with a 
numerical algorithm called the Parameter Switching algorithm which represents an elegant and easy way to numerically 
approximate any attractor of a dynamical system, belonging to a general class of systems, starting from a set of accessible 
parameter values which are switched in relative short period of times. The switching rule can be modeled by some piecewise 



continuous function. The algorithm is useful for example in practice, when a desired parameter value cannot be directly 
accessed. Also, it can help to understand what happens in some real systems when the control parameter is switched by 
natural or imposed causes. 

The Parameter Switching algorithm differs from the known control/anticontrol algorithms, where the parameter is 
generally slightly modified following some very precise rules in order to modify the behavior of some trajectory. Our 
algorithm allows the choice of any deterministic or even random switching rule within a set of parameter values, the result 
being an attractor which belongs to the set of all existing attractors of the underlying system. 

The paper is organized on two main parts concerning theoretical aspects and applications respectively, as follows: 
Section [2] presents the attractors synthesis, where the Parameter Switching is detailed, Section [3] presents the numerically 
evidence of the Parameter Switching algorithm convergence, and finally, in Section [4] four representative examples are 
analyzed. Additional information are presented in five Appendices. 



2 Attractors synthesis 

In this section, the used notions, results, assumptions and the underlying Initial Value Problem (IVP) modeling a general 
class of systems, continuous or discontinuous with respect to the state variable and of integer or fractional order, are 
presented. 

2.1 Preliminaries notions and notations 

All the considered systems can be modeled by the following IVP 

^ f{x{t))+pBxit) + Cs{x{t)), x{0)^xo, tel (1) 

where p is a real parameter, q stands for the derivative order (for g = 1, we have the known standard derivative, while 
for g ^ 1 we have the so called fractional derivative: cfl/dt'^), f : M" M" is a nonlinear vector valued function, 
at least continuous with respect to the state variable, / = [0, T] , T > , 5, C real n x n squared matrices, and 
s : M" M", s{x) — (si(a;i), . . . , s„(xAr))* is a piece-wise continuous function, being composed in most general cases by 
signum functions: Si{xi) = sgn{xi), i = 1, 2, . . . , 71, or e.g. step (Heaviside) functions. 

It is classically assumed that q G (0, 1]. Function of q and C, we can have the situations presented in Table [l] 
Throughout this paper the following assumption will be considered 

(HI) The IVP |Ip admits a unique solution (e.g. Lipschitz continuous). 

The control parameter p is considered to be a (periodic) piecewise constant function p : / — M (an example for a 
periodic function p is presented in Figjl]). As we shall see next, p can be a periodic or non-periodic function. Other form 
of functions for p can be found in ^TI] . 

Due to the piece-wise continuity of p, the IVP ([T]) becomes non-autonomous. However, for the sake of simplicity, next 
the time variable t will be omitted unless necessary. Therefore, the IVP can be written as follows 

d'^x 

— = fix) + pBx + Cs{x), x(0) = a;o, t e /. (2) 

Remark 1. The existence and uniqueness conditions for IVPs modeling DI and DF systems differ from those for CI 
systems, and are not presented here (for our class of DI systems they can be found in e.g. p[2j, while for differential 
equations of fractional order in 13J or [14 .) 

The systems chosen to represent in this work the four classes in Table [T] are three-dimensional, but the algorithm and 
the underlying results are applicable for any finite lower or higher dimension n. 

The examples treated here are presented in Table [2j Lorenz system, Sprott system [TS] , Lii system lTF] and a fractional 
variant of the Chua's system [T7] . Other examples can be found in [TJ [31 151 IHl [7] . 
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A global attractor, roughly speaking, is viewed in this paper as a state space region of a dynamical system where the 
system can enter but not leave, and containing no smaller regions (see e.g. [H]). The global attractor contains all the 
dynamics evolving from all initial conditions. In other words, it contains all the solutions, including stationary solutions, 
periodic solutions, as well as chaotic ones. The term of local attractor is used sometimes for attractors which are not 
global attractors [TOj . 

The global attractors may contain several local attractors. Therefore, a global attractor can be considered as being 
"composed" of the set of all local attractors for a given parameter p value and initial conditions. Each local attractor 
attracts trajectories from a subset (basin of attraction) of initial conditions (for details on the notions of local and global 
attractors we refer e.g. to [20 l [2T | [22 l l23]). 

Remark 2. (i) For the sake of simplicity, when a global attractor is composed by several local attractors, only a single 
local attractor will be considered (the choice can be made by appropriate selections of the initial conditions). Therefore, 
hereafter, by attractor one understands, simply, either one of the local attractors or the single local attractor which 
composes the global attractor; 

(ii) Due to the predominant numerical characteristics of the present work, without a significant loss of generality, the 
attractors will be considered as approximations, after neglecting a sufficiently long period of transients i24j . of the uj — limit 
sets (the set of points that can be limit of subtrajectories). Despite the fact that usually these sets are uncomputable, 
they can be numerically approximated. Therefore, in this paper the attractors are considered as being the plots of the 
Lj — limit sets. 

Let us use throughout the review the following notations 

Notation 1 - A the set containing the attractors depending on p, including attractive stable fixed points, limit cycles 
and chaotic attractors; 

- V the set of all p admissible values; 

- V^f = {pi,p2, ■ ■ ■ ,Pn} cV a. finite ordered subset of V; 

- An = {^1, • ■ • J ^at} C a the set of the attractors corresponding to V^j-; 

N 

- I — y ( y lij), where the adjoint subintervals lij are of time length m^h, where the "weights" are some positive 

] = 1,2.... 1=1 

integers, h > 0, for i = 1,2,...,N and all j (see Figj2] for the particular case of the first set of time-intervals In for 
^ = 1,2,3,4); 

- p* the average parameter 

N 

J2 "l^iPi 



1=1 

- A* the average attractor, obtained for p ~ p* . 

Remark 3. Taking account to the Assumption HI it follows naturally to define a monotone bijection F : Vn ~^ An for 
some fixed N. Therefore, to each p e Vj^/ corresponds a unique element A e Aj\f and reversely, for each A G Aj\f there 
exists p^Vm such that A = F{p) (Figj|]). 

2.2 Parameter Switching algorithm 

To prove that any attractor can be approximated by switching the parameter while the underlying IVP is integrated, we 
need a numerical algorithm to implement the switches that we name Parameter Switching (PS) algorithm. 
Let fix for some iV, the set Vn- Then, p* given by (|3|, can be rewritten in the following form 

N IN 

P* = ^ aiPi with ai^mi / ^ "i^, pi £ Pn- (4) 
1=1 I 1=1 
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N 

Because < 1 and at ~ 1 , p* enjoys the following property 

i=l 

PI. For every set Vn , P* given by ^ is a convex combination of pi, i = 1,2, . . . , N . 

To implement the PS algorithm, we have to integrate the IVP ^ with a numerical scheme for ODEs with single 
step-size h (e.g. the standard Runge-Kutta method). 

Let first consider that p is a periodic function of period Tq, i.e. p{t + Tg) = p{t) for all t in /. While the solution to 
the IVP ^ is numerically approximated, the parameter p is periodically switched within Vn in every consecutive time 
interval 7^-, following some designed scheme, denoted hereafter by Sh 

Sh = [pi P2 |/2j, • ■ • ,PN ] , j = 1,2, . . . , (5) 



which means that while the IVP ([2| is integrated, in each interval lij, p will be replaced by pi for every j = 1,2,.... Thus, 
for t G III , pit) = pi, for t S I21, p{t) = P2 and so on until /7V1, when p(t) = p^- On the next interval I12, again p{t) = pi 
and so on until the interval In2, when p{t) — pn- The algorithm repeats on the next set of intervals Ii3,i — 1,2, ... ,N 

N 

and so on periodically, until t > T. In other words, p is a piecewise constant and periodic function of period Tq — h mi 
having the following expression (See Fig |l|2[ ) 

p{t)=p„fortehj, i = l,...,N, j = l,2,... (6) 

The length of the time intervals lij will be taken as multiple of h: length{Iij) = ruih for each j. Therefore, for a fixed h, 
Sh can be noted in a simplified form 

Sh = [mipi , m2P2 , • • • , ttinPn] , (7) 
which means the following p infinite sequence 

TOiPi, TO2P2, • ■ • , ruNPN, Wipi, m2P2, • • ■ , m^PN, ■ ■ ■ 

For example, Sh = [2pi,p2] for a given h, means that for the time-interval of length 2h, p — pi then for the next time- 
interval of length h, p — p2 . Next, for two integration steps, p — pi, then for one integration step, p — P2 and so 
on 

(i.e. periodically with period Tq = {mi + "^2)^ = 3ft.). Schematically, Sh can be written as the infinite sequence 

[2pi,P2] = Pl,Pl,P2,Pl,Pl,P2, ■ ■ ■ 

The pseudocode of the PS algorithm is given in Table [Sj 
It is easy to verify that 

t+To 

P* = ^ J P{T)dT, t £ I. 
t 

Notation 2 Let denote by A° the attractor, obtained with the PS algorithm, called hereafter the synthesized attractor. 

Remark 4. It is easy to see that, for some given p, the relation ([s]) considered as equation, may have several solutions. 
For example, if we set N = 2, and want to obtain p* — 4 using the scheme Sh = [mipi, m2P2], for pi = 2 and p2 = 6, we 
can choose rui = 7712 = 1 but also mi = m2 = 3 to verify ([s]). If we fix mi = 3 and m2 = 1, in order to obtain p* = 4, we 
can use pi = 2 and p2 = 10, but also pi = P2 = 4. 



4 



3 Numerical proof of PS algorithm convergence 



In this section we prove numerically that for a chosen set of attractors An, the synthesized attractor A° obtained with 
the PS algorithm belongs to and, moreover, A° is approximatively identical to A* . 
In order to compare two attractors, we have to provide the following criterion 

Criterion We shall say that two attractors are approximatively identical (AT) if their trajectories in the phase space 
approximatively coincide, and the Hausdorff distance (Appendix [B|) is small enough. 

In our numerical experiments Hausdorff distance was of order of 10^'* — 10^'^. 

Due to the bijectivity of F, considering the total order over the set Vn, it is reasonable to consider that the following 
property holds 

P2. An is an ordered set endowed with the Vn order induced by the function F. 

Moreover, the same order can be found over the sets Vn and considered as intervals: [pi,pn] and [Ai, An] respectively 
and, without losing generality, we can consider that A^ — F{pi), i = 1,2, . . . , N (FigQ. This property is outlined in all 
bifurcation diagrams. 

Notation 3 Let denote AI by " " . 

Next, in order to prepare the proof of the main result regarding the parameter switching, we introduce the following 
lemma 

Lemma 1. Given N and Vu, A° = A* . 

Proof. The lemma has been checked numerically with tools such as: histograms, Poincare sections, time series, cross- 
correlation (Appendix [A| and Hausdorff distance (Appendix [b|) . The numerous examples, show that the attractor A° , 
obtained with PS algorithm and A* obtained for p ~ p* , are AI, the degree of the identity depending less or more on the 
system characteristics and, unavoidably, on the numerical errors. Hausdorff distance, for all considered systems, was of 
order of 10"'' - lO^^. □ 

The sketch of the analytical proof of this lemma, presented in [11] for the case of CI systems, can be found in Appendix 

o 

Remark 5. Applying the symbolic computation for several examples of CI systems, with the scheme Sh for < 3, the 
IVP ([2]) was integrated with the forward Euler method. The result shown that the (Euclidean) difference between the two 
solutions corresponding io p = p* and to p switched with the PS algorithm, is of order of 0{h^), the same as the error of 
the considered Euler method. 

Yet, the main result which can be numerically proved, can be introduced. 
Theorem 1. Civen N and Vjs/, A° belongs to {Ai,An)- 

Proof. By the properties PI and P2, it follows that A* G [Ai, An)- Next, by the Lemma[l] the attractor A° synthesized 
with some scheme Sh, is AI to A* . Thus, A° ^ A* and therefore A° belongs to {Ai, An), which completes the proof (see 
Figjs]). □ 

Summarizing, for every finite set Vj\f and numbers m^, the synthesized attractor A° will belong to {Ai,An), and differs 
from every attractor Ai S Aj\r, i — 1,2,. . .N (due to the convexity property). Reversely, any attractor of Aj^r can be 
considered as being synthesized with the PS algorithm, by means of a finite set of attractors of Ajs/. 

For the continuous case, the analytical proof in [TT], shows that the solutions of the equation ([2| with p switched within 
Vj\f with PS algorithm and that with p replaced with p* can be arbitrarily close. Therefore, the underlying invariant sets 
(attractors in our case) are also arbitrarily [AI) close ([E], Ch. 6). 
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Due to the mentioned convexity property, whatever kind of combinations of pi and values are considered for a fixed 
N , p* will still belong to the interval {pi,pn)- Therefore, for p we can chose any kind of piecewise continuous functions, 
with the only condition that their values and p* belongs within Vj^ (see examples in pTj). 

The proof of the convergence (analytically or numerically verified) does not depends on the periodicity of p but only 
on the convexity of p* . Therefore, it is obvious that not only periodic schemes ([t]) can be used, but even random ones [T] . 
One of the simplest way to implement randomly the PS algorithm, once N is fixed, is to chose first pi and in some 



random manner, after which the PS algorithm is started (see the example of Sprott system, Subsection 4.3 1. Obviously, 
there are several other random ways such as: choosing randomly rrii and pi while the PS is running, or switching the 
order of pi and so on. Now, the averaged p* has to be determined with the following relation 

N 

E m'iP^ 

P = ^ : (8) 

E ml 

i=l 

where denote the number of times when pi is chosen by the algorithm for t G /. 

Remark 6. (i) The "structural stability" of the PS algorithm presents some obvious limitations due firstly to his numerical 
approach (some details and other related aspects can be found in [T]). For example, for relative large values of rm, the 
trajectory of A° could present some "corners" (a maximum difference between rrii should generally be about (20 25)h). 
The values for pi could be taken over the entire set 7-V without distinguishable differences between A° and A* . Excessive 
number of decimals for p* could lead too to some differences between the two attractors A° and A* . Even for large values 
for N, A° and A* still remain close each other; 

(ii) The cross-correlation and time series show an interesting characteristic: the trajectories corresponding to A° and A*, 
even in the phase space and time representations are AI, they are dephased in time (see cross-correlation in the figures); 

(iii) For chaotic attractors, the AI is obtained only "asymptotically" since the necessary time to fully approximate the 
attractor is, theoretically, infinite. 

The PS algorithm can be used to "control" or "anticontrol" dynamical systems modeled by the IVP ^ when some targeted 
parameter value cannot be accessed directly (see ^ ) . For this purpose, we have to choose pi , rrii and some scheme Sh to 
obtain the targeted value p* . However, while almost all known control/anticontrol algorithms "force" some trajectory to 
change its characteristics and behavior, the PS algorithm allows to obtain any desired already existing attractor of Aj\f. 



4 Applications 

This section is devoted to the applications of the PS algorithm to the four classes of dynamical systems (see Tables [l] and 
[2]) to synthesize attractors. In this purpose we have to choose N, Vj^f and Sh for each system, such that a desired value p* 
(which can be taken e.g. from the bifurcation diagram) is obtained. 

To apply the PS algorithm for CI systems, we used the standard Runge-Kutta method (with the step size h of order 
between 10~^ and 10~^, depending on the characteristics of the considered system), while for the discontinuous and 
fractional systems, we have chosen special numerical methods. The bifurcation diagrams, time series, histograms and 
cross-correlations were determined and plotted superimposed for the first state variable xi. The Poincare sections have 
been determined for the plane = const. Some bifurcation diagrams, like the one for the Sprott system and especially 
for the fractional Lii and Chua systems, require an extremely long computer time (see Appendix [d| ) . For discontinuous 
systems (of integer and fractional order), some 'corners' can be remarked, typical for these kind of systems (the solution 
for the underlying IVP are generally not smooth [12 [5S]). As stated before, the AI was verified for all the considered 
systems via superimposed phase portraits, Poincare sections, histograms, time series and also with cross-correlation and 
Hausdorff distance. For all the considered cases, the results lead to the same conclusion: Lemma [T] applies to all considered 
classes of systems. 
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4.1 Continuous dynamical systems of integer order 

The PS algorithm was tested on several examples of CI systems such as: Lorenz, Chen, Rossler, Rabinovitch-Fabrikant 
[T], Hindmarsh-Rose neuronal system 0, networks and Lotka-Volterra [7 . Here, we consider the representative case 
of the Lorenz system. 

For this class of systems, the IVP ^ has to be considered for the particular case q — 1 and C — Onxn, namely (Table [2]) 

X = f{x) +pBx, x{0) = xo, te I. (9) 
The used numerical method is the standard Rungc-Kutta with integration step-size h = 0.01. 

• Let first consider the scheme ([t]) for N = 2: [mipi, m2P2\ with pi — 90, p2 = 96, and toi — m2 — 1. Then p* , given 
by the relation ( [3|, i s p* = (toi xpi+m2 'x P2) / {1^1+1^2) = 93. Applying the PS' algorithm, the synthesized attractor 
A° (red plot, Fig|6j a) is a stable limit cycle that is A I with the average attractor A* (superimposed blue plot over 
A°) for p* = 93. The AI is emphasized in addition to the over-plot in the phase space, by the superimposed Poincare 
sections with the plane X3 = 130 (Figjoja) and superimposed histograms for the first state variable xi too (Figj6]: ). 
The cross-correlation (Fig|6]d) shows that the time series corresponding to A° and A* are AI, but dephased. This 
time-difference between the corresponding trajectories is better remarked from the time series corresponding to the 
first component xi (see Figjoje). 

• Let next consider the case N = 5 with the scheme [2pi, 3p2, 2p3, 4p4, 3^5] for pi = 125, p2 = 130, = 140, p4 = 144 
and P5 — 220. In order to facilitate the use of the PS algorithm, the bifurcation diagram will be used (Figj?]). Now, 
p* — 154, and the synthesized attractor A° is again a stable limit cycle which is AI with A* (Figjsjf) even Ai_4 are 
chaotic and only A^ is a stable limit cycle (Figjsja-e). Both attractors A° and A* are AI (see Fig[8]g-h wherefrom 
the AI property can be remarked). The time series being dephased (Figj8]i,j), the trajectories of the attractors A° 
and A* are AI, but time dephased. 

• If for the same scheme [2pi, 3p2, 2p3, 4p4, Sps] we choose ps = 166 instead ps — 220, the synthesized attractor A° is 
chaotic and AI with A* for p* — 142.428 (Figj9|. However, now the ^/ is only an almost identity (see Remark [6] 
(iii)). The Poincare section (Figj9]c) was obtained with the plane X3 = 145. The cross-correlation shows that the 
underlying trajectories of A° and A* are dephased. Because the trajectories are chaotic, the time series to underline 
this time-difference is irrelevant in this case. 

For all analyzed examples, the Hausdorff distance was of order 10~^. 
Other examples of CI systems can be found in (llj . 

4.2 Continuous dynamical systems of fractional order 

Fractional mathematical concepts allow to describe certain real objects more accurately than the classical "integer" 
methods. Examples of such real objects that can be elegantly described with the help of fractional derivatives displaying 
fractional-order dynamics, may be found in many fields of science and exhibit a wide range of rich dynamics. Therefore, 
the fractional calculus starts to attract increasing attention of mathematicians but also of physicists and engineers (see 
e.g. EillZllllIll]). 

Many CF systems, can be modeled by the IVP ^ with q e (0, 1) and C = Onxn- The fractional derivative ^ is 
generally denoted using the Caputo differential operator of order q, Dl (see e.g. [30])- Thus, the IVP ([2| becomes 

Dlx^ f{x)+pBx, x('=)(0)==4'=\ (fc = 0,l,...,M -1). (10) 

[.] denotes the ceiling function that rounds up to the next integer, and 13™ = -^p^ , with m = \q] , is the standard 
differential operator of the integer order \q] E N. The Caputo operator, with starting point 0, has the following expression 

r[m - q) Jo 
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where F is the Gamma function (Appendix [d|). Because I?* has an m-dimensional kernel, m initial conditions need to be 
specified. Therefore, for the common case chosen in this paper < g < 1, we have to specify just one condition, in the 
classical form [31]: x^^\Q) — xq. 

To implement the PS algorithm in this case, it is necessary to choose a numerical method for the solution to the IVP 



(10). In this purpose we use the fractional Adams- Bashforth- Moulton method (see Appendix [D]) introduced in [5Tj . 
Let choose for our purpose the fractional variant of the Lii system (see Table[2]) which unifies the Lorenz and Chen systems, 
presented by Lii et al. in [35] . As many of the real fractional systems have the order of the fractional differential operators 
less than 1, we fix in this paper q = 0.9 (see [33]) which is a typical value exhibiting all the relevant phenomena (the 
dynamics of this system, as q varies, can be found in [16], while some aspects of the attractors synthesis of the fractional 
Lii system is presented in |34j). 



Now, the IVP ( 10 ) was integrated with the fractional Adams-Bashforth-Moulton method with step size h = 0.005 and 



15000 20000 steps, in function of the dynamics of the synthesized attractor A° 



A chaotic attractor A° can be obtained with the scheme [Ipi, \p2\ (see Fig 10 c) with pi = 32 and p2 = 34.5. The 



attractors corresponding to pi and p2 are plotted in Fig |10| a, b which, as can be seen in the bifurcation diagram 
in Fig 11 are stable limit cycles, p* — 33.25, and due to the chaotic behavior, the AI between A° and A* is only 



asymptotic (see the histograms in Fig 10 d and the Poincare sections in Fig 10 f). Even A° and A* are chaotic, from 



the cross-correlation (Fig 10 e), one can see that they are dephased, but still AI. 



If we choose the scheme [lpi,lp2], with pi — 33.5 and p2 = 35.5, a stable limit cycle A° , for which p* — 34.5, 
is obtained (Fig 12 1. The AI can be remarked from the Poincare section and histograms (Fig 12 b,c). Again, the 
time difference between the underlying time series can be seen from the cross-correlations (Fig 12 d) and time series 
(Figfl2]e) 



For all tested CF systems, the Hausdorff distance was only of order of 10~^, compared e.g. with CI systems, where it was 
of order of 10""^. This is explainable due to the well known 0{h?) error bound for the used one-step Adams-Bashforth- 
Moulton method for fractional systems (detailed discussions on errors can be found in 13 ). 

4.3 Discontinuous dynamical systems of integer order 

Differential equations with discontinuous right-hand side, model a whole variety of realistic applications: dry friction, 
electrical circuits, oscillations in visco-elasticity, brake processes with locking phase, oscillating systems with viscous 
dumping, electro-plasticity, convex optimization, control synthesis of uncertain systems and so on (see e.g. [3S1 [Ml I3Z] 
and the references therein). 

For our class of DI systems, (? = 1 and C ^ Onxn and the IVP ([2| becomes 



X — f{x) + pBx + Cs{x), x{0) — Xq, t£l. 
In this case, the right-hand side is discontinuous for a null set of points M where s vanishes, and continuous in Z? = . 



Obviously, the IVP (11) cannot be solved with classical methods. For example, for the equation 

X — 2 — 3sgn(a::), 



where M — R\{Di IJ D2) = {0} with Di = (—00, 0), D2 = (0, 00), the classical solutions, for x ^ 0, are 

xit)-- 



5t + Ci, X <0, 
-t + C2, x > 0, 



(11) 
(12) 



(13) 



with the integration constants Ci, C2 but, as t increases, these solutions tend to the line x — 0, where they cannot continue 
to evolve along this line since the function x{t) = does not satisfy the equation (Fig[l3|). Thus, there is no classical 
solution starting from 0. 



^In [38j a classification of systems modeled by the IVP is presented. 
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Therefore, the problem has to be restarted as a differential inclusion by using, for example, the Filippov regularization 
(Appendix [e|). Thus, the IVP (111 transforms into a differential inclusion (set-valued IVP) 



X £ f{x) +pBx + CS{x), x(0) = xq, for almost all t € I, 
where S is the setvalued variant of s. On mild assumptions, a differential inclusion has a solution that happens to be 



even unique, but it could have multiple solutions too. To find them numerically, in our particular case of the IVP (11), 
we can use the standard Runge-Kutta method within D and a special numerical method for differential inclusions in M 
(the simplest forward Euler method here, see Appendix [e]). 



Once we set the numerical method for the IVP (11), we can apply the PS algorithm. For this class of DI systems, we 
choose the Sprott system [15 (Table [2]). 

• First, let us choose N — 2 and the scheme [lpi,lp2] for Pi = 0.5 and p2 = 0.528, for which the corresponding 



attractors Ai and A2 are chaotic (Fig 14). We have chosen this scheme such that the obtained average value 



p* = 0.514 belongs to a stable periodic window in the bifurcation diagram. Therefore, the synthesized attractor 



is a stable limit cycle (Fig 15 c). The AI is underlined by the superimposed Poincare sections (with the plane 



X3 — 0, Fig 15 d) and histograms (Fig 15 e). The time-difference between the trajectories is remarked from the 



cross-correlation (Figjl5|f) and time series (Fig|l5|g). 

As seen in Section 3, PS algorithm can be implement in random manners too. For example, for N = 100 if one 
choose randomly rrii G {1,2,3} and pi G [0.45,0.65] for i = 1,2,..., 100 with uniform distribution, and with the 



obtained values we launch PS, the synthesized attractor A° is still AI to the average attractor A* (Fig 16 1. However, 
taking into account the asymptotic generation of chaotic attractors, and the relative large value for N, the small 
differences between the two attractors, A° and A* are explainable. 

4.4 Discontinuous dynamical systems of fractional order 

There are real discontinuous dynamical systems which display fractional-order dynamics. We consider here the following 
class of DF systems, modeled by the IVP ([2| for C 7^ Onxn and p < 1 

cPx 

— = f{x) + pBx + Cs{x), x(0) =xo,t£l. (14) 
In |39j it is proven that the IVP admits solutions which can be numerically determined. 



Shortly, the IVP is transformed first into a differential inclusion via the Filippov regularization (as in Subsection 4.3). 
Next, using the Cellina's theorem (see e.g. [40l p. 84]) the set-valued IVP of fractional-order is transformed into a 
continuous single-valued of fractional-order IVP (see for continuous approximation of DI systems the way chosen in 
|41j). The approximation is made in a sufficiently small £- neighborhood of the discontinuity points. To be precise, let us 
consider the simplest example of the scalar function, widely used in examples: s{x) — csgn{x). To approximate s(x) in 
an ^-neighborhood oi x = 0, we can choose one of the simplest function, the sigmoid hg 



For our general case of the IVP ( 14 ) , the continuous approximation leads to the following continuous IVP of fractional-order 



d^x ( Cs{x), for x(^M, , 

— -/(x)-pi?.= | ^^(^^^ /or XGM, (16) 

where h^[x) is the e-approximation of Cs{x) in the e-neighborhood of the points x € M, verifying the continuity condition 
h^{x) — Cs{x) on the boundary of the £-neighborhood [H]. In this way, the discontinuous IVP became a continuous one of 
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fractional order and a numerical scheme for fractional-order differential equations, such as the Adams-Bashforth-Moulton 



method presented in Subsection 4.2 can be used. 

We consider for this case the fractional variant of a discontinuous Chua system presented by [17] (Table [2]) for q = 0.98 
(smaller values gives not rich dynamics). 



As it can be seen from the bifurcation diagram (Fig 17 1, for p e (12, 12.55), there exists a narrow band of a kind of 



"chaotic saddle" . Within this window, the underlying chaotic attractors look as being "embedded" within this transient 



chaos (see for example the attractor Ai in Fig 18)) 



• Using the scheme [2pi,p2] with pi = 12.5 and p2 = 17, the obtained synthesized attractor A° is AI with A* for 
p* = 14 (Fig[l8|. 

As shown in [33] . we found numerically that in these systems of lower than third-order (i.e. 3*9 which, for q < 1, is less 
than 3) chaos still may appear (as it is known, in the case of integer order, according to the well-known Poincare-Bendixon 
theorem, chaos appears only at systems of minimum order three). 



Conclusions 

In this review we have presented the parameter switching algorithm according to which any attractor of a dynamical 
system belonging to a large class of systems, may be numerically approximated (synthesized). The attractors synthesis is 
achieved by using the PS algorithm, which switches periodically or randomly the parameter. This facility is enabled by 
the convexity of p. The average and synthesized attractors are AI and their underlying trajectories, time dephased. The 
review is legitimated by the more than ten published papers each of them containing several applications. 

As expected, the performance of the PS algorithm is limited due to the errors of the used numerical method, the 
length of the time-subintervals Ik, k = 1,2, ...,N, the number of digits for p, the step size h and the distance in the 
parameter space between pk- Thus, we found that TV is not a critical parameter (it could be even about 100). The length 
of Ik (i.e. the value for m^) is a critical parameter indicating for how long time the control parameter of the considered 
system can take the values p — p^. We found that a maximum value for can be taken about 25/i. For p, 3 — 4 digits 
are enough to be compatible with the smallest distance between the pk in the bifurcation diagrams. Moreover, some 
real physical chaotic systems may have an infinite number of different states or limit cycles with infinite period. But a 
computer simulated system has a finite number of states; if the precision of the computer is n bits and the system to be 
modeled has k variables, the total number of system states is limited to 2*^*"; hence, given a determined state, it will be 
repeated sooner or later and the system will become periodic, with a period equal to the separation of the two occurrences 
of the state. The PS method can alleviate this inaccuracy and make possible the approximation of a computer simulated 
system to a real one, although it may be necessary to use a sequence of parameter values lasting as the whole segment of 
the system to be modeled (se e.g. [H]). 

Some open problems are: the analytical proofs for the Lemma [T] for CF, DI and DP systems, not only for CI as done 
in [TT]; an analytical proof for the continuity of the bijection F; a detailed study of the time delay between the trajectories 
of A° and A*\ the effect of noise on the results; a comparison with the complex systems (fractals), where the parameter 
switching may lead also to some interesting results [S]. 



APPENDIX 



A Cross Correlation 

As known, the cross-correlation of two signals is a measure of the similarity of two waveforms. The cross-correlation has 
ranges from -1.0 to -1-1.0. The closer it is to -1-1 or —1, the more closely the two compared variables are related. The 
correlation of two signals (the attractors underlying trajectories in our case) may indicate that one of them is delayed in 
time with respect to the other. The maximum value (close to unity) of this cross-correlation is obtained when the two 
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signals are in closest alignment with each other. The value —1 means the signals are identically matched but opposite in 
phase, while a value approaching zero indicates a low degree of similarity (see the blue band around the horizontal axe 
in our images). In this paper, the results were obtained with the crosscorr Matlab function with approximate 95 percent 
confidence interval. 



B HausdorfF distance 

The HausdorfF distance in a metric space, measures how far two compact nonempty subsets are from each other. The 
classical Hausdorff distance between two (finite) sets of points, A and B, is defined as [43, p. 114] 



Dh{A, B) — max < sup inf d{x, y), sup inf d{x, y) 
where d(x, y) is the classical distance between two points in the considered space. 

If the two sets are curves, Dh is defined as the maximum distance to the closest point between the curves. Thus, if the 
curves are defined as the sets of ordered pair of coordinates A = {oi, 02, . . . , a„} and B = {61, &2, • ■ • , bm} the distance to 
the closest point between a point to the set B is 

d{ai, B) — min \\bj — ai\\ . 
3 

Thus, the HausdorfF distance is 

di/(A, i?) — max < max {d(ai, -B)} , max \d{hj, A)} 

C Sketch of the analytical proof of Lemma [l] 

Next, the main steps of the proof presented in pQ] for the lemma, ensuring the AI between A° and A* in the case of CI 
systems, are pointed out. 

Consider the IVP ^ with C = 0„xn and q = I satisfying the assumptions stated in Section [2] and expressed for the 
general case of M" , in the following form 

i{t) ^ f{x{t))+p{t/X)Bx{t), a;(0) = xo,i e / = [0,00), (17) 

where A G is a positive real number which will be stated later, and p : I ^ M" is considered as a piecewise continuous 
periodic function with period To, and mean value q, i.e. 



1 



t+To 

p{u)du ~ q, t G I. 



Let us define the average model of (17 1 obtained with the PS algorithm, expressed as follows 

y = f{y) + qBy, y(0) = yo- (18) 



The IVP (17) models the PS algorithm and generates the synthesized attractor A° , while the IVP (18) represents the 
system whose solution approximates the average attractor A*. 

We have to prove that the solutions of the equations (fl7|) and ( 18 ) differ by less than for A sufficiently small, via the 



so called order function defined in terms of approximation^ 
Let next suppose that (18) satisfies the assumption HI and admits s : / ^ M as the unique solution. 



^The order function i5(A^), introduced in |44l p. 11], implies that there exists k s.t. |5(A'^)| < kX^ when A is sufficiently small. 
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Linearizing (18) on a neighboriiood of s, one obtains the following IVP 



e{t) = [E{t) + qB]e{t), e(0) = eo, 
where e{t) = y{t) — s{t) and E{t) denotes the Jacobian of / evaluated at s{t). 



(19) 



Because s{t) is the solution in (18), e{t) = for i S /, should be a solution of (19). 
If we linearize the IVP ( 17) for a: G Tg (the domain of attraction of e = 0) one obtains 



3{t)^[E{t)+p{t/X)B]e{t), e(0) = 



Co, 



where e{t) = x{t) — s{t). 

Then, the theorem ensuring the AI between the attractors of the dynamical system modeled by the IVP (17 1 and IVP 
( 18 ) can be enounced. 



Theorem Let assume that Eq. ( 19 ) is uniformly exponentially stable, i.e. there exist the constants if > 0, /i > such 
that 

e{t) < X||£o||exp(-^i). 

Then, for eo — Eq, there exists a positive scalar A > 0, such that hmt_i.oo \ \e{t) — e{t)\\ — (5(A^), where 6{X^) is an order 
function. 

Proof. The complete proof can be found in [TT| and it mainly follows the idea given in Chapter 4 of ^44J- The existence 
interval / is partitioned as follows: / = [0, XT] \J [XTq, 2XTo] ■ ■ ■ . In each subinterval = [uATq, (n + l)ATo] , n = 1,2,..., 
Eq. (19) has the solution e„(f). If on these subintervals, the initial condition is chosen e„(nATo) = e(nAro), using a 
generalized Peano-Baker series |45], the Gronwall's inequality, through straightforward algebraic operations, the following 
inequality is inductive proved 

||e((n + l)ATo) - e„((n + l)ATo)|| < ^(A^), 



for any n. Taking the limit n — )■ oo, the proof is complete. 



□ 



D Adams-Bashforth-Moulton scheme for fractional ODEs 



Next, a brief presentation of the Adams-Bashforth-Moulton scheme [3T] is presented. Let consider the IVP (10). Specifi- 
cally, the method implies a discretization of I with grid points ti = hi, i = 0, 1, . . . with a preassigned step size h. First, 
a preliminary approximation a;f_^x ^'^^ ^{ti) (the predictor phase) is computed via the formula 

j=o i=0 

where have the form 

Then, the final approximation Xi+i for a;(ti+i) (the corrector phase) is 

r-ji-i hi f ' 

with 

ao,^+l^^^+^ -{i-q){i + l)\ 
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and 



for j = l,2,...,i. 

The Gamma function, F, is approximated in this work with the so-called Lanczos approximation |46j 

r(^) = ^LoZ!£l(;^ + 5.5)-+o.5g-(-+5.5)^ 

for z € C with Re(z) > 0. The coefficients pi are shown in Table |4j 

While in the standard methods for ODEs of integer order, the current approximation Xk depends only on the results 
of a few backward steps, like all reasonable numerical methods for fractional differential equations, the fractional scheme 
Adams-Bashforth-Moulton requires the entire backward integration history at each point in time. Thus, each calculated 
value Xk depends on all previous values xo,xi, . . . , Xk-i- This characteristic implies a serious drawback with respect to the 
required computing time. For example, to obtain 4000 points within some attractor, about 8 x 10® iterations are necessary. 
However, this is necessary to appropriately reflect the memory effects possessed by fractional differential operators. 
A detailed analysis of this method can be found in |13j and a background on fractional differential equations is presented 



E Filippov regular izat ion 

Let consider the following general IVP with discontinuous right-hand side 

i = /(x), x{0)=xo, tel, (20) 

with / locally bounded on M". In particular, the discontinuity is due to the discontinuity of the state variable, of the 
associated vector field, of Jacobian (partial derivatives) or higher order discontinuity. The continuity domain consists in 
a finite m number of open regions Di C K", i — 1,2, . . . ,m, the discontinuity set M being M = K"\ [J^-^ Di. 



The IVP (20) may have not any solutions in the classical sense. Therefore, in this paper we have chosen the way given by 



137] I which transforms the IVP (20 1, via the differential inclusion approach, into a multi- valued Cauchy IVP 

X £ F{x), x(0) = Xq, for almost all t £ I, (21) 

where, F : M" =| M" is a set-valued function into the set of all subsets of M". For our class of systems, F can be defined 
using the so called Filippov regularization 

F{x) — con lim f{x'), 

where con means the convex hull and lima;'->.a; f{x') represents the set of all limits for all convergent sequences f{xk) with 
Xk X. For X e M, F{x) is a set, while for x ^ M, F{x) consists in a single point f{x). As example, for the sign function 
the Filippov regularization gives the following set-valued function 

{-1}, a:<0, 
Sgnix) = { [-1,1] a: = 0, 
{+1} X > 0. 



For example, the Filippov regularization applied to the Example ([12| leads to the following differential inclusion 

ie 2-3Sgn(a;). (22) 



Definition [A7\ A generalized solution (or Filippov solutions) of the IVP (20) is an absolutely continuous vector- valued 



function x : / — >■ i?" verifying the IVP (21 ) for a.a. t £ I 
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Even the IVP (20) may have no classical solutions, the setvalued IVP (21) may have a unique or several generalized 
solutions. For example, the equation (12), after regularization becomes the set-valued problem (22) and has the following 

-t + xq for t < xq and x(t) = for t > xq- In other words, the solution can 



0. If Xq < 0, the solution is x(t) ^ 5t + x^ for t < x'q and x(t) = ior t > Xq 



generalized solutions: if xo > 0, then x{t) 
be prolonged continuous along the axis x 
(Fig{l3]) 

The background of differential inclusions and their solutions can be found e.g. in |40j and [48] . The existence and 
uniqueness of solutions for our class of DI systems are presented in |49j and will be not considered here. 

To solve numerically the IVP (20) special numerical methods for differential inclusions are necessary. However, for our 
class of IVPs, due to the presence of s functions, the discontinuity appears only in a finite null set M, where actually the 
IVP is a set-valued problem. For the points x G Di, the IVP is a continuous problem. Therefore, we can integrate in Di 
the IVP (20) using e.g. the standard Runge-Kutta method, while for x S M, a numerical method for the corresponding 
differential inclusion x € F(x) has to be used. Precisely, when the trajectory enters the discontinuity surface, we have 
to choose for derivative of solution, generally for some finite time, a value within the set F{x) while a numerical method 
for differential inclusions is used (the simplest one is the adapted forward Euler method see e.g. I25j ). For example, 
when a; = in the Example (12), we have to solve the differential inclusion x G [—1,5]. There are several possibilities to 
manage this problem using e.g. so called selection strategies (see [SO]). In this paper we used the simplest way, namely 
the random strategy which implies a randomly choice of a value within F{x) (the interval [—1,5] in our example). There 
are several possibilities to find the moments when the trajectory enters and leaves the discontinuity surfaces (see e.g. [S7] ) 
during which the chosen method solves the differential inclusion. 
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Tabic 1: Classification of the considered dynamical systems modeled by the IVP 



Continuous of Integer order (CI) systems 
Continuous of Fractional order (CF) systems 



Discontinuous of Integer order (DI) systems 
Discontinuous of Fractional order (DF) systems 



Table 2: Systems utilized in this paper. 



Type Order System 
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Table 3: Pseudo-code of the PS algorithm 



CHOOSE Sh, T, h 
REPEAT 

FOR i = 1 to 

FOR fc = 1 to 

one step integration of the IVP ([2j) for p = pi 
t = t + h 
ENDFOR 
ENDFOR 
UNTIL t > T 
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Table 4: CoefRcients of the Lanczos approximation. 



i Pi 

75122.6331530 

1 80916.6278952 

2 36308.2951477 

3 8687.2452971 

4 1168.9264948 

5 83.8676043 

6 2.5066283 
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Figure 1: Piecewise constant periodic function p : / — ^ M (sketch). 
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Figure 4: Order induced by F in A 
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Figure 5: Attractor synthesis: sketch of the proof of Theorem [TJ 
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Figure 6: Scheme [mipi,m2P2] with mi = 90, m2 = 96, and pi = P2 = 1 appHed to the Lorenz system, p* = 93. (a) A° 
and A*; (b) Poincare sections; (c) Histograms; (d) Cross-correlations; (e) Time series. 
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Figure 7: Bifurcation diagram for the Lorenz system. 
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Figure 8: Scheme [2pi, 3p2, 2p3, 4p4, Sps] for pi — 125, p2 = 130, = 140, = 144 and p^ = 220 applied to the Lorenz 
system, p* = 154. (a)-(e) The attractors Ai corresponding to p^, i = 1, . . . , 5; (f) A° and A*; (g) Histograms; (h) Poincare 
sections; (i) Cross-correlations; (j) Time Series. 
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Figure 9: Same scheme as in Figjs] [2pi, 3p2, 2p3, 4p4, Sps] but with ps = 166 instead P5 = 220. p* = 142.428. (a) A° and 
A*; (b) Histograms; (c) Poincare sections; (d) Cross-correlations. 
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Figure 10: Scheme [Ipi, lp2] with pi ~ 32 and P2 — 34.5 appUed to the fractional Lii system, p* = 33.25. (a),(b) The 
attractors, Ai and A2; (c) The attractors A° and A*; (d) Histograms; (e) Cross-correlations; (f) Poincare sections. 
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Figure 12: Scheme [lpi,lp2], with pi — 33.5 and p2 = 35.5 for the Lii system, p* — 34.5. (a) A° and A*; (b) Poincare 
sections; (c) Histograms; (d) Cross-correlations; (e) Time series. 
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Figure 14: Bifurcation diagram for the Sprott system. 
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Figure 15: Scheme [Ipi, IP2] for pi — 0.5 and p2 — 0.528 applied to the Sprott system, p* — 0.514. (a),(b) Ai and A2; (c) 
A° and A*; (d) Poincare sections; (e) Histograms; (f) Cross-correlations; (g) Time series. 
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Figure 17: Bifurcation diagram for the fraetional Chua system. 
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Figure 18: Scheme [2pi,p2\ with pi = 12.5 and p2 — 17, apphed to the fractional Chua system, p* = 14. (a),(b) Ai and 
A2; (c) A° and A*; (d) Poincare sections; (e) Histograms; (f) Cross-correlations. 
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